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EXECUTIVE  SUMMARY 


SRI  International  modeled  the  comminution  and  flow  of  ceramic  under  conditions  like 
those  in  ceramic  armor  at  the  nose  of  an  advancing  penetrator.  We  have  delivered  a  model, 
FRAGBED2,  that  models  fracture,  comminution,  compaction,  and  fragment  flow,  and  can  be 
implemented  in  hydrocodes.  Fragment  motion  is  treated  by  analogy  with  atomic  dislocation 
theory.  That  is,  blocky  fragments  are  imagined  to  glide  in  small  increments  along  fixed  planes  in 
the  material.  As  fragments  flow  one  way,  lines  of  holes  flow  the  other  way.  Compaction  occurs 
by  the  outflow  of  holes,  and  dilatancy  occurs  by  the  influx  of  holes,  and  dilatancy  occurs  by  the 
influx  of  holes.  Comminution  is  treated  by  an  “overstress”  type  of  rate  model.  The  rate  of 
comminution  is  zero  below  a  threshold  stress,  saturates  above  a  saturation  stress,  and  varies 
smoothly  in  between. 

FRAGBED2  is  an  improvement  over  models  used  presently  for  armor  simulations 
because  it  captures  the  relevant  physics.  FRAGBED2  parameters  are  directly  related  to  readily 
measured  material  properties,  such  as  the  fracture  toughness,  friction  coefficient,  relative  size  of 
pre-existing  flaws,  density,  and  porosity. 

We  tested  FRAGBED2  by  comparing  results  from  three  types  of  experiments:  spherical 
cavity  expansion  tests,  thick-walled  cylinder  collapse  tests,  and  ballistic  tests.  In  the  cavity 
expansion  tests,  an  explosive  charge  in  a  spherical  cavity  machined  in  a  split  block  of  ceramic 
produced  a  distribution  of  fragment  sizes  as  a  function  of  distance.  By  examining  the 
distribution  posttest,  we  obtained  parameters  for  calibrating  the  comminution  model.  The  thick- 
walled  cylinder  test,  conceived  by  Nesterenko,  was  modified  by  introducing  a  taper  that  allowed 
a  complete  range  of  strains  to  be  produced.  Extensive  shear  banding  was  observed  in  recovered 
specimens.  Two-dimensional  hydrocode  simulations  with  FRAGBED2  were  able  to  model  the 
final  deformation  of  the  cylinder,  but  were  not  able  to  model  the  observed  shear  banding.  We 
found  that  three-dimensional  flow  in  the  specimen  makes  it  difficult  to  compare  these 
experimental  results  with  the  results  of  two-dimensional  calculations.  Thus,  three-dimensional 
hydrocode  calculations  probably  will  be  required  to  test  the  ability  of  FRAGBED2  to  model 
shear  banding.  FRAGBED2  did  a  good  job  of  matching  the  ballistic  test  results.  The 
experimental  results  were  for  long-rod  penetrations  into  alumina  ceramic,  and  were  obtained 
from  the  literature. 
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SECTION  1 

INTRODUCTION  AND  BACKGROUND 


SRI  International  modeled  the  comminution  and  flow  of  ceramic  under  conditions  like 
those  in  ceramic  armor  at  the  nose  of  an  advancing  penetrator.  Penetrating  weapons  or  explosive 
charges  in  brittle  materials  such  as  ceramics,  concrete,  or  hard  rock,  produce  fracture  and 
fragmentation  near  the  cavity  boundary  to  form  a  bed  of  fragmented  or  granulated  material. 
Penetration  occurs  by  further  comminution  of  the  material  into  a  finely  granulated  bed  (called  the 
Mescall  zone)  and  the  subsequent  flow  of  the  granules  out  of  the  way  of  the  advancing  penetrator 
[1],  as  shown  schematically  in  Figure  1.  In  quasi-static  and  dynamic  tests  under  conditions  of 
compression  and  shear,  the  yielding  behavior  of  brittle  frictional  materials  can  be  interpreted  as 
the  result  of  frictional  sliding  of  debonded  granules.  Thus,  modeling  the  formation  and  flow  of  a 
granulated  bed  is  key  to  computing  cratering  and  penetration  in  brittle  materials,  and  to 
interpreting  laboratory  tests  of  such  materials. 


Figure  1 .  Damage  pattern  in  ceramic  during  long-rod  penetration. 
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In  a  prior  paper  [2],  we  presented  a  granulated  material  model  for  use  in  finite  element 
hydrocodes  applied  to  penetration  of  ceramic  armors  or  hard  rock.  The  model,  named 
FRAGBED,  included  a  mesomechanical  description  of  shear  flow  and  associated  dilatancy.  It  is 
a  non-local,  multiplane,  plasticity  model  (see,  for  example,  Batdorf  and  Budianski  [3], 

Curran  et  al.  [4],  and  Bazant  et  al.  [5])  that  proved  useful  in  computational  simulations  and 
associated  interpretations  of  penetration  experiments  in  which  ceramic  armors  were  attacked  by 
long-rod  penetrators  [1]. 

However,  a  drawback  to  FRAGBED  was  that  the  pore  compaction  was  handled  from  a 
continuum  viewpoint,  whereas  the  shear  flow  and  associated  dilatancy  were  treated  from  a 
mesomechanical  viewpoint.  In  fact,  both  dilatancy  and  pore  compaction  arise  from  the  same 
mesomechanical  processes.  In  addition,  the  important  process  of  granule  comminution  was  not 
treated  in  FRAGBED,  and  the  fracture  processes  that  initially  produce  the  fragmented  bed  were 
treated  only  cursorily. 

We  expanded  the  FRAGBED  model  to  be  a  complete  micromodel  for  the  fracture, 
fragmentation,  comminution,  shear  flow,  dilatancy,  and  pore  compaction  processes  in  ceramics. 
We  call  this  expanded  model  FRAGBED2,  and  the  original  model  FRAGBED  1. 


2 


SECTION  2 
APPROACH 


The  FRAGBED  models  are  mesomechanical  and  average  the  behavior  inside  a  relevant 
volume  element  (RVE),  as  discussed  by  Curran  et  al.  [4],  Bazant  [5],  Nemat-Nasser  [6],  and 
numerous  other  authors.  The  RVE  must  contain  many  granules.  As  described  in  [2],  the 
FRAGBED  approach  to  modeling  the  flow  of  granulated  material  is  to  focus  not  on  the  granules 
themselves,  but  rather  on  the  holes  between  the  granules.  Figure  2  provides  a  schematic  view  of 
this  approach. 

By  analogy  to  atomic  dislocation  theory,  a  hole  large  enough  to  allow  sliding  of  a  granule 
into  it  is  called  a  vacancy,  and  strings  of  such  holes  are  called  dislocations.  Just  as  for  the  atomic 
case,  the  dislocations  can  be  mobile  or  pinned,  will  have  edge  and  screw  components,  and  can 
glide  or  climb.  The  dislocations  can  become  pinned  when  obstacles  stop  their  motion,  and  can 
later  become  unpinned  by  granule  rearrangement.  As  discussed  in  [2],  this  analogy  is  useful 
because  it  allows  us  to  easily  relate  non-elastic  slip  in  the  granular  bed  to  macroscopic  plastic 
strain  rate.  By  casting  the  model  in  the  framework  of  multiplane  plasticity  theory,  the  analogy 
between  granular  flow  and  slip  in  single  crystals  is  direct,  and  many  of  the  techniques  and  results 
of  the  community  engaged  in  mesomechanical  modeling  of  crystalline  plasticity  can  be  applied. 
We  can  thus  use  terms  commonly  associated  with  atomic  dislocation  theory  to  refer  to  granular 
flow  processes. 

Before  the  dislocation  analogy  can  be  applied,  the  material  must  be  converted  from  an 
intact  material  to  a  fragmented  bed.  FRAGBED  1  currently  uses  a  very  simple  damage  evolution 
algorithm  based  on  a  combination  of  tensile  and  shear  strains.  When  the  damage  function  attains 
a  critical  value,  the  material  cohesion  and  tensile  strength  reach  zero  and  the  material  element  is 
declared  fragmented.  Thereafter  the  fragmented  bed  model  is  used,  with  the  average  fragment 
size  and  dislocation  density  specified  as  input  constants. 

In  FRAGBED2  we  introduced  three  changes.  First,  FRAGBED2  provides  a  simple 
model  of  the  coalescence  of  cracks  nucleated  at  flaws  to  form  an  initial  granule  size  distribution 
(FRAGBED  1  assumed  a  constant  initial  granule  size;  this  option  remains  available  in 
FRAGBED2).  However,  this  initial  fragmentation  algorithm  remains  very  approximate.  The 
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Dislocation  motion 


Figure  2.  FRAGBED2  dislocation  analogy. 

(a)  Case  I  -  unimpeded  motion,  (b)  Case  II  -  impeded  motion, 
(c)  Definition  of  parameters,  (continued) 
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second  change  was  to  introduce  evolution  equations  for  the  average  fragment  size  and  the 
dislocation  density,  i.e.,  we  introduced  a  comminution  model.  The  third  change  was  to  introduce 
a  micromechanical  model  for  porosity  evolution.  This  model  replaced  the  dilatancy  model  and 
the  continuum  model  for  compaction  used  in  FRAGBED1.  Both  the  dilatancy  and  pore 
compaction  are  a  consequence  of  dislocation  flux  across  the  boundary  of  the  RVE;  dilatancy  is 
caused  by  dislocations  flowing  into  the  RVE,  whereas  compaction  is  caused  by  dislocations 
flowing  out  of  the  RVE. 
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SECTION  3 

DESCRIPTION  OF  FRAGBED2 


FRACTURE  AND  FRAGMENTATION  PROCESS 

We  assume  that  initially  intact  brittle  material  is  converted  to  a  fragmented  bed  by  the 
activation,  growth,  and  intersection  of  cracks  that  initiate  at  inherent  flaws  in  the  material.  For 
example,  in  an  initially  intact  ceramic  under  compressive  and  shear  stresses,  one  could  imagine 
that  the  flaws,  which  may  consist  of  relatively  weak  grain  boundaries,  microcracks,  or  other 
stress  raisers,  will  become  active  when  the  shear  stresses  exceed  the  cohesion,  or  when  local 
anisotropy  and  inhomogeneity  cause  local  tensile  stresses  even  under  global  compression. 

[Under  uniaxial  strain  conditions,  this  would  occur  at  the  Hugoniot  elastic  limit  (HEL)].  In  a 
jointed  hard  rock  formation,  the  same  process  could  occur  under  compressive  and  shear  loads 
that  vary  slowly  over  the  scale  of  the  joint  spacing.  For  common  ceramics  or  intact  hard  rock, 
this  will  require  compressive  stresses  on  the  order  of  10  GPa,  but  would  occur  at  much  lower 
pressures  for  heavily  jointed,  weakly  cemented  rock.  For  the  ceramic  penetration  and  tapered 
cylinder  experiments  to  be  discussed  later,  we  assume  that  the  primary  fragmentation  process 
near  the  nose  of  the  penetrator  and  in  the  tapered  cylinder  is  grain  boundary  debonding.  Thus,  the 
initial  fragment  size  distribution  is  that  of  the  ceramic  grains.  (In  these  cases,  significant  shear 
flow  occurs,  and  it  is  the  fragment  sizes  of  the  flowing  material  that  are  relevant.) 

For  cavity  expansion  experiments,  although  the  above  compression  and  shear  process  will 
occur  in  the  immediate  vicinity  of  the  cavity  boundary,  the  fragmentation  process  farther  from 
the  cavity  boundary  occurs  via  the  intersection  of  tensile  cracks.  These  tensile  cracks  are 
nucleated  by  the  spherically  divergent  shock  wave,  which  produces  tensions  normal  to  the 
direction  of  propagation.  These  larger  fragments  exhibit  negligible  shear  flow.  We  next  develop 
a  simple  algorithm  for  this  fragmentation  process. 

A  complete  description  of  this  process  would  begin  by  using  a  mesomechanical  model 
such  as  that  of  Simons  et  al.  for  concrete  [7],  to  calculate  the  fragment  pattern.  For  a  recent 
review  of  such  approaches,  see  Rajendran  [8].  However,  we  postpone  such  complexity  by  using 
a  simple  algorithm  as  a  preprocessor  to  establish  the  initial  fragmented  bed,  as  discussed  below. 
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In  cases  of  penetration  or  cratering,  we  observe  experimentally  that  cone  cracks  radiate 
into  the  target,  and  these  are  intercepted  by  spoke  cracks.  The  spoke  cracks  are  sketched  in 
Figure  3.  The  shaded  region  is  the  loaded  region,  or  cavity.  Both  spoke  and  cone  cracks  arise 
from  membrane  tensions  produced  by  the  spherically  divergent  shock  wave  and  associated 
particle  velocity  field.  Surface  observations  commonly  show  many  such  cracks  near  the  load 
center  and  fewer  cracks  at  relatively  large  distances  from  the  center. 


Figure  3.  Radial  spoke  crack  pattern  on  surface  of  brittle  target. 

The  compressive  stress  in  the  shock  propagation  direction  is  followed  by  tesions  that  can 
cause  lateral  cracks  that  intersect  the  spoke  and  cone  cracks  to  produce  a  fragmented  bed.  Also, 
the  spoke  and  cone  cracks  branch  and  meander  to  intersect  and  form  fragments. 

The  divergent  shock  wave  will  produce  high  hoop  tensions  at  high  strain  rates  (>104  s'1) 
in  the  target  material  near  the  load  center.  These  tensions  will  initiate  many  spoke  and  cone 
cracks,  and  their  average  spacing  can  be  estimated  with  the  Grady  algorithm  (for  a  recent  review, 
see  [9])  for  brittle  material: 
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s 


pc 


fd£v' 


(3-D 


where  S  is  the  crack  spacing,  Kc  is  the  fracture  toughness,  p  is  the  density,  ds/dt  is  the 
continuum  strain  rate  at  the  time  of  fracture,  and  c  is  the  longitudinal  sound  speed. 

The  strain  rate  of  interest  is  the  rate  at  which  the  material  is  stretching  in  the  membrane 
direction  in  the  shock  pulse,  which  is  given  by 

(3-2 

dt  r 

where  u  is  the  radial  particle  velocity  and  r  is  the  distance  from  the  impact  site.  We  obtain  u 
from  the  jump  condition 


P 

u  =  — 
pc 


(3-3) 


where  P  is  the  compressive  shock  stress  and  c  is  the  shock  velocity.  Combining  Eqs.  (1)  -  (3) 
yields 


The  shock  pressure  in  an  elastic  fluid  decays  according  to  the  equation 


P  = 


(3-4) 


(3-5) 
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where  a  is  the  initial  cavity  radius,  which  we  can  approximate  with  the  penetrator  radius,  and  Pq 
is  the  impact  pressure.  Combining  Eqs.  (4)  and  (5),  and  setting  S  equal  to  the  initial  fragment 
size.  Bo,  yields 


Bn 


5 Kc 

\P°a  J 


.X 


(3-6) 


As  discussed  later,  in  modeling  the  cavity  expansion  test,  we  used  the  above  equation  to  create 
the  initial  fragmented  bed  before  we  began  the  hydrocode  simulations.  We  adjusted  the 
comminution  parameters  to  match  experimental  observations.  Thus,  although  Eqs.  (3-6)  give 
initial  sizes  somewhat  larger  than  the  grain  sizes,  the  comminution  process  soon  reduces  the  sizes 
to  below  the  average  grain  size. 


COMMINUTION  PROCESS 

We  return  to  the  simplified  picture  of  the  fragmented  bed  discussed  in  [2],  as  shown  in 
Figure  2.  We  assume  an  initial  size  distribution  of  fragments,  idealized  here  as  equisized  square 
blocks. 

A  key  result  of  confinement  at  boundaries  is  that  dislocation  flux  across  material  element 
boundaries  is  decreased,  thereby  inhibiting  non-elastic  flow  of  the  fragmented  bed.  Perfect 
confinement  would  in  fact  prevent  long-rod  penetration  of  a  ceramic  target  because,  although  the 
fragmented  bed  could  flow  slightly  into  any  initial  porosity,  the  granules  would  ultimately  have 
no  place  to  go. 

We  next  consider  the  comminution  of  the  initial  fragment  bed  of  Figure  2.  We  assume 
that  the  comminution  process  mainly  consists  of  blocks  pushing  on  each  other  to  produce  local 
shear  and  tensile  stresses  that  in  turn  cause  the  blocks  to  fracture. 

Figure  2  shows  that  fracture  of  a  single  block  to  form  two  smaller  blocks  causes  the  small 
holes  above  and  to  the  right  of  the  fractured  block  to  become  a  single  dislocation  (a  dislocation  is 
a  line  of  holes  big  enough  to  allow  non-elastic  block  motion),  because  as  the  top  part  of  the 
fractured  block  moves  it  produces  a  hole  twice  as  big  to  its  left,  thereby  allowing  the  full-sized 
blocks  to  move.  The  lower  part  of  the  fractured  block  moves  to  the  left  at  the  same  time  to  form 
a  right-moving  dislocation.  Thus,  the  fracture  of  the  block  nucleates  a  pair  of  edge  dislocations 
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traveling  in  opposite  directions,  and  the  dislocation  nucleation  rate  is  directly  tied  to  the 
comminution  rate. 

If  all  the  Case  II  blocks  were  to  break  in  two  in  the  same  manner,  the  situation  would 
become  that  of  Case  I,  with  no  block  interference,  but  with  twice  the  number  of  dislocations.  In 
all  cases,  the  nucleation  process  as  such  does  not  cause  increased  porosity  (the  hole  area  to 
volume  ratio  remains  the  same).  A  generalization  of  this  process  will  be  discussed  later. 

DILATANCY  AND  PORE  COMPACTION  PROCESSES 

The  internal  dislocations  in  Figure  2  move  up  and  left  and  down  and  right  to  cross  the 
RVE  boundaries,  and  result  in  the  hole  closures  and  associated  compaction.  At  the  same  time, 
dislocations  from  adjacent  RVEs  may  move  into  the  RVE,  resulting  in  dilatancy.  Thus,  the 
dilatancy  or  pore  compaction  is  a  result  of  the  net  flux  of  dislocations  across  RVE  boundaries. 
This  process  will  be  discussed  in  more  detail  in  Section  4. 
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SECTION  4 

EVOLUTION  EQUATIONS 


NON-ELASTIC  SLIP 

As  discussed  in  [2],  the  Orowon  [10]  equation  is  used  on  each  slip  plane.  That  is,  on  the 
ith  slip  plane, 


M. 

dt 


(4-1) 


where  }  ?  is  the  non-elastic  (plastic)  slip  strain  on  the  ith  plane,  Ndi  is  the  mobile  dislocation 
density  on  that  plane,  2?,  is  the  block  (granule)  size  in  the  direction  of  slip,  biBi  is  the  size  of  the 
dislocation  hole  in  the  direction  of  slip  (a  sort  of  macro  Burger’s  vector),  vdi  is  the  dislocation 
velocity,  vhj  is  the  corresponding  block  velocity,  and  g{  is  a  variable  coefficient  that  depends  on 
granule  geometry.  A  detailed  discussion  is  given  in  [2]. 

As  in  [2],  the  total  strain  rate  on  a  slip  plane  is  the  sum  of  the  elastic  and  non-elastic  rates: 


dt 


Jl(21l)+2Zl. 
2g|^  dt  J  dt 


(4-2a) 


dt 


1 

dip 

K( l-f») 

^  dt  j  dt 

(4-2b) 


where  the  porosity  <p  =  0(.  (see  later  discussion),  ji  is  the  total  shear  strain  on  the  ith  plane, 

l 

1  (.  is  the  resolved  macroscopic  shear  stress  across  the  plane,  £ v  is  the  total  volumetric  strain,  P  is 
the  pressure  in  the  porous  material,  <pt  is  the  porosity  (the  ratio  of  void  volume  to  total  volume) 

associated  with  the  ith  plane,  and  G  and  K  are  the  shear  and  bulk  moduli  of  the  intact  material, 
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respectively.  The  division  by  (l  -  <p)  in  the  first  term  on  the  right-hand  side  (rhs)  of  (2b)  is  based 
on  the  assumption  of  a  partial  pressure  relationship: 


p=(\-<f)p. 


(4-2c) 


where  Ps  is  the  compressive  mean  stress  in  the  intact  material.  We  have  also  neglected  terms  of 
the  order  of  PjKi) -(f).  In  the  applications  to  be  discussed,  the  porosity  is  restricted  to  a  few 
percent  by  comminution  and  associated  compaction,  and  this  approximation  is  justified.  Our 
model  is  not  suitable  in  any  case  for  large  porosity,  because  in  that  event  the  dislocation  approach 
appears  unrealistic. 

POROSITY 

We  assume  for  simplicity  that  all  porosity  is  associated  with  dislocations.  The  total 
porosity  (p  is  defined  as  the  sum  of  the  porosities  (pt  associated  with  each  plane,  where  (p-t  is 
defined  as  the  ratio  of  void  volume  associated  with  slip  on  the  i  th  plane  to  the  total  volume.  To 
simplify  the  notation,  the  i  subscript  will  be  dropped  in  further  equations,  and  the  equations  will 
be  assumed  to  refer  to  the  i  th  plane  unless  otherwise  specified.  Thus,  the  value  of  porosity  on 
each  slip  plane  associated  with  Nd  mobile  dislocations  per  unit  area  is 

.  (p  =  nNdbB 2  (4-3) 


where  we  have  approximated  the  individual  mobile  dislocation  area  as  nbB2 ,  where  n  in  Figure 
2  is  unity.  For  more  realistic  fragmented  bed  geometries  than  those  of  Figure  2,  we  can  expect  n 
to  be  greater  than  unity  to  allow  unimpeded  block  motion.  We  have  arbitrarily  chosen  n  =  3/2 . 

The  porosity  can  change  only  by  hole  migration  across  a  material  element  boundary 
(dislocation  flux).  That  is,  as  discussed  in  [2],  the  non-local  nature  of  FRAGBED  arises  from  the 
choice  of  a  control  volume  (RVE)  that  contains  many  granules  and  associated  dislocations. 
Applying  the  divergence  theorem  to  the  dislocation  density  for  this  control  volume  yields  (on  a 
single  slip  plane) 
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(4-4) 


dNj_  ajf/.vj) 

dt  dt  dx 


where  the  first  term  on  the  right  refers  to  the  net  nucleation  rate  of  mobile  dislocations  (the 
nucleation  rate  minus  the  pinning  rate),  and  the  second  term  on  the  right  refers  to  the  net  flux  of 
dislocations  across  the  surface  of  the  material  element.  As  seen  in  Figure  2,  the  nucleation  and 
pinning  processes  do  not  change  the  porosity.  Thus,  only  the  second  term  on  the  right  of 
Eq.  (4-5)  is  used  to  obtain  the  evolution  relation  for  the  porosity,  and  for  constant  b  and  B 


_  j2  d(A TjVj.  ) 

2  dx 


(4-5) 


where  x  is  in  the  direction  of  slip  on  the  i  th  plane. 

The  parameter  b  is  important  because  it  helps  determine  the  porosity  through  Eq.  (4-5). 
Its  value  is  determined  physically  by  the  granule  geometry  and  configurations.  In  the  current 
form  of  FRAGBED2,  however,  the  value  of  b  is  simply  taken  to  be  an  adjustable  constant  with 
magnitude  between  zero  and  unity.  (A  value  larger  than  unity  would  mean  that  the  dislocation 
jog  would  be  greater  than  the  granule  size.) 

GRANULE  AND  DISLOCATION  VELOCITIES 

We  also  need  an  equation  for  the  dislocation  velocity,  which  is  in  turn  related  to  the 
granule  velocity  [see  Eq.  (4-1)].  We  can  obtain  the  form  of  the  equation  from  the  simplified 
picture  of  Figure  2.  Consider  a  cubic  block  of  size  B  undergoing  shear  stress  r  across  an 
interface  that  may  produce  a  slip  in  the  x  direction.  The  frictional  resistance  to  slip  is  JLion.  The 
impulse-momentum  relation  is 


(' t-/X(Tn)B2dt  =  psB3dv 


(4-6a) 


where  //  is  the  coefficient  of  intergranular  friction,  on  is  the  compressive  normal  stress  across 
the  interface  between  granules,  ps  is  the  solid  density,  and  v  =  dx/dt  is  the  block  velocity. 
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Integration  of  Eq.  (4-6a)  gives 


v  = 


f2(r-//<7„> 

'  PsB 


(4-6b) 


From  Eqs.  (4-6a)  and  (4-6b),  the  time  for  the  granule  to  move  across  the  gap  bB  is 


t  =  B 


2  Pi 

t-pcn 


(4-6c) 


When  the  block  moves  a  distance  of  bB,  the  dislocation  moves  a  distance  B.  Thus,  the 
dislocation  velocity  vd  is  given  by 


T-//CT 

2  p,b 


n 


0<vd  <c 


(4-6d) 


In  [2],  a  more  approximate  derivation  was  used  that  produced  a  linear  dependence  on 
(z  -pon).  The  new  formulation  Eq.  (4-6d)  produces  weaker  dependence  of  dislocation  velocity 
on  stress  and  hence  lower  strain  rate  sensitivity  than  the  FRAGBED1  formula. 


DISLOCATION  NUCLEATION  PROCESS 

As  discussed  above  and  illustrated  schematically  in  Figure  2,  fracture  of  a  granule  has  the 
effect  of  nucleating  new  mobile  dislocations  (two  of  opposite  sign  in  the  example  of  Figure  2). 

In  an  actual  fragmented  bed,  the  granules  will  not  be  equisized,  but  will  have  a  distribution  of 
sizes.  We  could  replace  the  equisized  blocks  in  Figures  1  and  2  with  a  commonly  observed 
Poisson  granule  size  distribution: 


Ntl(R)=Nb, 


-RIB 

e 


(4-7) 
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where  Nbg  (R)  is  the  number  of  granules  per  unit  volume  with  radii  greater  than  R,  Nbl  is  the 
total  number  of  granules  per  unit  volume,  and  B  is  thus  the  characteristic  size  of  the  granule  size 
distribution. 

We  could  then  reinterpret  the  Bs  in  the  previous  equations  as  applying  to  the  B  in 
Eq.  (4-7).  However,  in  the  results  presented  in  this  report,  we  have  continued  to  assume  that  the 
granules  in  a  computational  cell  all  have  one  size,  B. 

The  net  nucleation  rate  of  mobile  dislocations  can  be  related  to  the  comminution  rate 
dB/dt  through  Eq.  (4-3).  Because  the  nucleation  process  does  not  create  porosity,  we  can 
differentiate  Eq.  (4-3)  with  respect  to  time,  and  set  dQjdt  =  0,  to  obtain 


9AT,  _  2N„  dB 
dt  B  dt 


(4-8) 


Substituting  Eq.  (4-8)  into  Eq.  (4-4)  yields 


3N„  2NtdB  B(NjVd)  (49. 

dt  B  dt  dx 

GRANULE  COMMINUTION  PROCESS 

We  imagine  a  situation  like  that  depicted  in  Figure  2.  Fragments,  or  rather  the  voids 
between  them,  can  be  thought  of  as  dislocations  gliding  along  planes  oriented  in  discrete 
directions.  Our  present  concern  is  the  further  fracture  of  the  fragments,  such  as  might  occur  at 
the  potential  fracture  site  indicated  in  Figure  2. 

In  the  comminution  model,  fragment  size,  B(x,t),  decreases  at  a  rate  determined  by  the 
amount  the  driving  stress  exceeds  a  threshold  defined  by  the  current  fragment  size  and  the 
fracture  toughness.  At  sufficiently  high  stresses,  the  rate  saturates.  We  imagine  that  fragments 
break  in  two  once  every  period  At.  Then 


3B  HB  2B  lBf 
ot  At  At  2  J 


(4-10) 
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where  /  =  l/A t  is  the  frequency  of  breaking.  We  assume  that  fragments  contain  crack-like  flaws 
of  size  tjB  and  that  linear  elastic  fracture  mechanics  applies.  Fragments  break  if  the  stress 
intensity  induced  by  the  driving  stresses  exceeds  the  fracture  toughness,  i.e.,  if 


Od  > 


(4-11) 


The  frequency  is  postulated  to  reach  a  limit  because  the  cracks  cannot  grow  faster  than  some 
speed,  CR ,  possibly  but  not  necessarily  the  Rayleigh  wave  speed.  The  limit  breaking  frequency 
is  then  CR/ 2 .  We  ignore  the  possibility  for  crack  branching.  Thus,  we  expect  the  comminution 
rate  to  behave  as  shown  in  Figure  4. 

We  can  model  this  behavior  using 


oB 

ot 


(4-12) 


where  S(pd,  B)  is  a  shelf  function  that  makes  a  transition  from  0  to  1  as  the  driving  stress 
parameter  Od  increases  from  below  the  critical  threshold  to  above  the  saturation  limit.  The 
threshold  is  determined  by  the  current  fragment  size,  B,  and  the  saturation  limit  is  chosen  to  be 
some  fixed  multiple  of  the  threshold. 

From  an  infinity  of  possibilities  we  have  chosen  the  cumulative  distribution  function  for 
the  Weibull  distribution  [1 1]  for  our  shelf  function.  Specifically, 


dB 

dt 


' 

“ 

ecr0 

1  „ 

1-exp 

(a  A 

2<tw 

=-  -c  A 

— 

a 

> 

2 

w 

(4-13) 


where  e  is  the  base  of  the  natural  logarithm,  the  parameter  ow  specifies  the  half-width  of  the 
transition  and  the  parameter  o0  anchors  the  location  of  the  transition.  See  Figure  4.  The  e  arises 
because  we  normalize  the  Weibull  shelf  function  so  that  the  slope  at  o0  is  2ow/\  (the  width 
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divided  by  the  rise  of  1).  Rather  than  specify  the  width  and  location  of  the  transition  directly, 
however,  we  choose  to  specify  the  threshold  stress  for  the  beginning  of  the  transition  and  a 
parameter  m  >  1  such  that  the  saturation  stress  is  m  times  the  threshold  stress.  The  threshold 
depends  on  the  current  fragment  size,  the  fracture  toughness,  and  the  relative  sizes  of 
microcracks.  The  threshold  stress,  as  given  by  the  right-hand  side  of  Eq.  (4-1 1),  is  simply  the 
result  obtained  from  linear  elastic  fracture  mechanics,  with  any  geometrical  factors  accounting 
for  fragment  and  crack  shape  and  orientation  suppressed. 

The  width  of  the  transition  is  defined  as  the  difference  between  the  threshold  and 
saturation  stresses.  The  location  of  the  transition  is  the  mean  of  the  threshold  and  saturation 
stresses.  Thus, 


and 


=  (* 


m 


\4tjB 


(4-14) 


<7o 


nK) 

4tjB 


\l/2 


(4-15) 


At  a  constant  driving  stress  and  with  an  initial  fragment  size  that  places  the  threshold  below  it, 
the  comminution  rate  oB/ot  will  decrease  toward  zero,  as  can  be  shown  by  substituting  Eqs.  (4- 
14)  and  (4-15)  into  the  differential  equation  (4-13)  and  recasting  in  non-dimensional  form.  In 
effect,  the  rate  curve  in  Figure  4  shifts  to  the  right  as  fragment  size  decreases. 

The  comminution  equation  in  the  non-dimensional  form  is 


dB_ 

ot 


< 

1-exp 

- 

f  e[m+l]  "N 

jg8(m-l] 

_ 

l  ) 

_L 

(4-16) 


where 


B  = 


16  VO]  B 
7rKf(m+if 


(4-17) 
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2i 


and 


gQ7gj 

xKf(m  +  lf 


(4-18) 


As  B  approaches  zero  from  above,  dBj cH  approaches  zero,  and  comminution  ceases.  For  the 
parameter  values  in  Table  1, 


B(o)=  1.10xl05m15(0) 


(4-19) 


and 

t  =  2.75xlOV‘t  (4-20) 

i.e.,  for  an  initial  fragment  size  around  10  pm,  the  non-dimensional  initial  size  is  around  1,  and 
for  a  duration  of  1  ps,  the  non-dimensional  duration  is  around  275. 


Table  1 

COMMINUTION  MODEL  PARAMETERS 


Parameter 

Symbol 

Value 

Relative  crack  size 

n 

0.1 

Maximum  crack  speed 

5000  m/s 

Driving  stress 

4  GPa 

Fracture  toughness 

4  MPa/m 

Saturation  stress  multiplier 

m 

1.15 

We  can  exercise  the  FRAGBED2  comminution  equations  independently  by  considering 
the  compression  and  shearing  of  a  thin  layer  of  ceramic  that  is  infinite  in  extent.  If  we  assume 
that  stresses  are  homogeneous  and  equilibrated  (i.e.,  that  the  layer  is  thin,  then  the  comminution 
equation  (4-16)  is  a  set  of  ordinary  differential  equations  in  time,  one  equation  per  slip  plane.  In 
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general,  the  comminution  law  differential  equation  (4-16)  cannot  be  integrated  analytically,  but 
we  can  solve  numerically  for  fixed  stress  using  Mathematica  software  [12]. 

Several  numerical  solutions  and  their  exponential  decay  approximations  are  compared  in 
Figure  5  for  a  small  matrix  of  the  parameters  5(0 )  shown  in  Table  2.  We  performed  L2D 
hydrocode  [13]  calculations  using  a  single  element  to  check  our  implementation  of  the 
comminution  model.  The  Mathematica  numerical  solution  and  the  L2D  solution  are  virtually 
identical. 


Table  2 

MATRIX  OF  COMMINUTION  MODEL  PARAMETER  VARIATIONS 


TYl  5(0)  — > 

0.5 

1.0 

1.5 

1.05 

X 

X 

X 

1.10 

X 

X 

X 

1.15 

X 

X 

X 

For  very  large  times,  the  solution  to  the  comminution  equations  under  constant  stress 
loading  approaches  a  power  law  decay,  independent  of  5(0) : 


8(m-l) 

j§(?)°c  (4-21) 

where  e  is  the  base  of  the  natural  logarithm.  For  certain  values  of  5(0) ,  the  solution  follows  a 
power  law  decay  for  all  times.  For  other  values  of  5(0) ,  the  solutions  approach  power  law  decay 
more  slowly  for  larger  values  of  m. 

LOCAL  STRESSES 

Local  stress  variations  around  the  average  or  continuum  stress  state  occur  within  a  RVE, 
and  these  local  stresses  are  important  in  the  comminution  process.  Following  an  approach  by 
Costin  [14],  local  stresses  are  calculated  on  each  plane  that  are  different  from  the  average  or 
continuum  stresses.  These  local  stresses  are  meant  to  reflect  the  existence  within  the  RVE  of 
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local  stress  variations  around  the  average  due  to  granule  interactions.  The  local  stresses  are  used 
to  govern  the  comminution  within  a  RVE,  but  are  not  used  to  transfer  momentum  to  adjacent 
meshes.  (The  average  stresses  are  used  in  the  equations  for  conservation  of  mass  and  energy.) 

In  the  comminution  equations,  we  therefore  allow  od  to  be  a  local  stress  defined  by 

a j  =  /tt,T2  +  m2<x2  (4-22) 

where  the  m’s  are  adjustable  parameters. 
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SECTION  5 

NUMERICAL  IMPLEMENTATION  IN  HYDROCODES 


SOLUTION  PROCEDURE  FOR  SLIP  PLANE  EQUATIONS 

We  follow  the  work  by  Peirce  et  al.  [15]  on  slip  in  single  crystals.  We  decompose  the 
total  deformation  tensor 


F=FeFp 


(5-1) 


into  an  elastic  deformation  tensor  and  a  plastic  deformation  tensor.  Given  the  plastic 
deformation  F^,  we  can  solve  for  the  elastic  deformation  Fe  from  Peirce’s  Eq.  (18).  The  total 
deformation  F  is  input  into  the  FRAGBED  routine.  Because  as  the  elastic  constitutive  equation 
is  assumed  to  be  isotropic,  we  can  make  Cauchy  stress  a  direct  function  of  \e  (see  Malvern  [16], 
Eq.  6.8.35),  where  Ye  is  defined  from  the  left  polar  decomposition  of  Fe: 

Fe=VeR,  (5-2) 


and  where  Rf  is  the  elastic  rotation  tensor  (see  [16],  Eq.  4.6.1).  We  use  logarithmic  strain 


£=Ln(Ve) 


(5-3) 


where  Ln  stands  for  the  tensor  logarithm  of  a  tensor.  Finally,  we  have  the  evolution  equation  for 
the  plastic  deformation  tensor: 


(5-4) 
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Here,  jk  is  the  non-elastic  shear  strain  rate  on  the  k  th  slip  plane,  and  sk  and  n*  are  the  slip 
direction  and  slip  plane  normal  unit  vectors,  respectively.  The  symbol  0  denotes  the  dyadic 
product.  Henceforth,  we  will  suppress  the  superscript  k  but  assume  that  it  is  there. 

Vectors  s  and  n  are  referred  to  as  the  plastic  configuration  (determined  by  Fp  ),  and  s*  and 
n*  are  the  corresponding  vectors  in  the  current  configuration  (determined  by  F). 

n*=nF;‘  (5-5a) 


The  slip  plane  normal  n  is  simply  a  constant  vector  given  by  the  original  orientation  of  the  plane. 
The  assumption  is  that  only  elastic  deformations  change  the  orientation  of  a  plane.  The  slip 
direction  s*  is  given  by  the  following  equations 


t  =  ncT 

(5-5b) 

an  =n  t 

(5-5c) 

s*  =  t  -  fr„n* 

(5-5d) 

s=f;V 

(5-5e) 

Here,  G  is  the  stress  tensor,  t  is  the  traction  vector  on  the  plane  and  <jn  is  the  normal  stress  on  the 
plane.  Equations  (5-5c,d)  ensure  that  vectors  s*  and  n*  are  orthogonal.  From  Eq.  (5-5a,e)  it 
follows  that  s  and  n  are  orthogonal  whenever  s*  and  n*  are  orthogonal.  Equation  (5-5e)  follows 
from  the  definition  of  a  deformation  tensor.  Equation  (5-5a)  does  not  have  a  similar  physical 
meaning,  but  is  simply  a  formula  that  generates  a  vector  that  is  normal  to  any  vector  that  lies  in 
the  slip  plane.  The  various  vectors  s,  s*,  n,  and  n*  must  be  normalized  before  use  in  equations 
like  (5-4). 

Analytically,  the  plastic  deformation  tensor  will  always  be  incompressible  for  the 
following  reason.  The  plastic  strain  rate  tensor  on  each  individual  plane,  }  kpsk  0  n* ,  is  traceless 
because  the  tensor  product  between  two  vectors  that  are  normal  to  each  other  is  traceless.  The 
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sum  of  a  set  of  traceless  tensors  is  also  traceless.  Equation  (5-4)  has  the  same  form  as  the  rate  of 
the  total  deformation  tensor,  see  [16],  Eq.  4.5.14: 


—  =  LF  (5-6) 

dt 

where  F  is  the  deformation  tensor  and  L  is  the  velocity  gradient.  Whether  we  choose  to  interpret 
the  sum  of  the  tensors  }ksk  0n*  as  a  plastic  velocity  gradient,  Lp,  or  not,  mathematically  the 
rate  of  the  plastic  deformation  tensor  has  the  same  form  as  the  rate  of  the  total  deformation 
tensor.  Therefore,  because  a  traceless  L  (zero  volumetric  strain  rate)  leads  to  incompressible 
flow,  a  traceless  Lp  will  lead  to  incompressible  plastic  flow. 

Numerically  there  is  potentially  a  problem  with  incompressibility.  We  use  the  following 
finite  difference  equation  to  update  Fp: 


(5-7a) 


Superscript  n  stands  for  time  step  level.  The  rhs  is  correctly  time  centered  at  n  + 1/2  because  Lp 
is  ordinarily  centered  on  half  time  steps.  Rearranging  we  get 


(5-7b) 


The  equation  can  be  solved  as  a  linear  equation  system  for  the  unknown  F"+1  using  Cramer’s 
rule.  Therefore,  we  again  rearrange 


j 


(5-7c) 
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Everything  on  the  rhs  is  known.  In  principle,  the  fact  that  the  determinant  of 
[l  -  (At/2)L;;  J  must  be  non-zero  poses  a  restriction  on  At.  In  practice,  a  zero  determinant  means 

that  the  flow  is  so  quick  that  the  cell  inverts  in  a  single  time  step.  In  that  case,  the  calculations 
would  abort  before  applying  the  above  equations. 

When  we  test  Eq.  (5-7c)  above,  we  find  that  incompressibility  is  satisfied  to  within 
machine  roundoff.  We  do  not  yet  completely  understand  how  this  fortuitous  situation  comes 
about,  but  the  correct  time  centering  of  the  rhs  of  Eq.  (5-7a)  may  be  critical. 

The  plastic  flow  is  modeled  entirely  as  slip  on  discrete  slip  planes.  We  have  equilibrium 
on  a  certain  plane  when 


%  <c  +  juon 


(5-8a) 


where  c  is  the  cohesion,  while  we  will  get  slip  when 

'i>c  +  [ion  (5-8b) 

When  we  have  slip  on  one  or  several  planes,  the  plastic  deformation  tensor  has  to  be  modified  to 
reduce  the  resolved  shear  stresses.  Because  the  dislocations  have  a  limited  mobility,  the 
reduction  in  the  resolved  shear  stress  is  subject  to  a  limit  on  the  plastic  strain  rate. 

7,  S  gbBNjV,  (5-9) 


where  g  is  a  proportionality  factor,  b  is  the  burgeris  vector,  B  is  the  block  size,  Nj  is  the 
dislocation  density,  and  is  the  dislocation  velocity.  An  initial  thought  is  to  enforce 

i=c  +  juon  (5-10) 


on  all  planes  that  want  to  slip. 
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The  following  example  with  a  pair  of  slip  planes  at  right  angles  shows  why  this  approach 
does  not  generally  work.  One  plane  is  along  the  x  axis,  the  other  plane  along  the  y  axis.  We 
assume  2D  symmetry,  which  means  that  the  direction  of  the  resolved  shear  stress  is  in  the  x-y 
plane.  On  the  plane  with  the  normal  in  the  y  direction,  Eq.  (5-10)  requires 

'iyx=c  +  fiO)y  (5-1  la) 


while  on  the  plane  with  the  normal  in  the  x  direction,  Eq.  (5-10)  requires 

txy=c+noxx 


(5-1  lb) 


The  symmetry  of  the  stress  tensor  requires  tyx  =  txy,  while  at  the  same  time  cr^  will  not 
generally  be  equal  to  <jyy.  This  example  shows  that  we  cannot  generally  expect  to  fulfill 
Eq.  (5-10)  on  several  slip  planes  at  the  same  time. 

We  therefore  implemented  a  scheme  where  we  solve  one  plane  at  a  time.  The  plane  with 
largest  plastic  strain  rate  is  processed  first.  The  processing  order  is  important  because  a  change 
in  shear  stress  on  one  plane  may  change  the  stress  state  on  some  other  plane.  We  must  also 
check  at  the  next  plane  to  be  processed  to  see  if  the  no  slip  condition,  Eq.  (5-8a),  has  been 
fulfilled.  If  not,  we  back-calculate  a  trial  strain  rate  from 


GAt 


(5-12) 


where  G  is  the  shear  modulus.  This  equation  gives  the  strain  rate  that  relieves  the  overshear  in  a 
single  time  step.  This  strain  rate  is  limited  by  Eq.  (5-9).  The  shear  modulus  enters  as  one  G 
because  we  are  slipping  on  one  plane.  The  cohesion  does  not  appear  in  Eq.  (5-12)  because  it  is 
set  to  zero  after  the  first  slip.  We  have  found  that  the  cohesion  makes  little  difference  to  the 
results,  so  we  have  not  implemented  anything  more  elaborate. 
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Equation  (5-12)  assumes  an  infinitesimal  linear  Hooke’s  law.  Because  of  constitutive 
and  geometric  non-linearities  we  have  to  iterate  while  fine  adjusting  y  .  The  ultimate  goal  is 


=0 


(5-13) 


All  this  is  subject  to  the  limit  on  the  strain  rate  from  Eq.  (5-9). 

The  model  assumes  that  there  is  a  porosity  <ps  associated  with  each  slip  plane,  given  by 

<ps=l.5NdbB2  (5- 14a) 

The  total  porosity  (f)  is  given  by  the  linear  sum  of  Eq.  (5-14)  over  all  the  slip  planes. 

In  the  present  formulation  the  porosity  affects  the  pressure  equation  of  state.  We  are 
using  a  Mie-Griineisen  P(V,  E)  relation.  The  specific  volume  V  that  enters  the  equation  of  state 
is  modified  to 


V  =  (l -<f>ys  (5- 15b) 

where  V  is  the  macroscopic  specific  volume.  At  zero  porosity  the  macroscopic  and  equation  of 
state  volumes  are  the  same.  As  the  porosity  increases  for  a  constant  macroscopic  volume  the 
pressure  increases,  too. 

WELL-POSEDNESS,  THERMODYNAMIC  RESPECTABILITY,  MESH  SIZE 
INDEPENDENCE,  AND  INTERNAL  CONSISTENCY  REQUIREMENTS 

FRAGBED2,  a  mesomechanical  internal  state  variable  model,  is  similar  to  the  type 
discussed  by  Coleman  and  Gurtin  [17],  who  provided  a  recipe  that  guarantees  well-posedness 
and  thermodynamic  respectability.  However,  FRAGBED2  does  not  exactly  fit  the  Coleman  and 
Gurtin  recipe  for  internal  state  variable  models  because  it  contains  a  flux  term  in  the  evolution 
equation  for  the  mobile  dislocation  density  Nd ■  However,  we  can  use  an  approach  developed  by 
Whitham  [18,  19]  to  analyze  the  shear  wave  propagation  equation  across  each  slip  plane,  and 
thereby  establish  well-posedness  and  thermodynamic  respectability. 

Consider  one  of  the  slip  planes  in  FRAGBED2.  If  we  define  the  slip  to  be  in  the  x 
direction,  and  the  normal  to  the  slip  plane  to  be  the  y  direction,  and  x  and  y  are  Lagrangian 
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coordinates,  then  the  conservation  of  momentum  and  the  mass  equations  on  that  slip  plane 
become: 


Momentum: 


du  _  dt 
P  dt  dy 


(5-16) 


Mass: 


d }  _  du 
dt  dy 


(5-17) 


where  u  is  the  particle  velocity  in  the  x  direction,  p  is  the  density,  7"  is  the  shear  stress  across  the 
slip  plane,  and  y\ s  the  strain  (displacement  gradient)  associated  with  the  slip. 

The  constitutive  relation  for  slip  on  the  plane  is 

(5-, 8) 

where  G  is  the  shear  modulus,  b  is  the  relative  macrodislocation  size  (dimensionless),  B  is  the 
granule  size,  Nd  is  the  mobile  macrodislocation  density  (number  per  unit  area),  and  vd  is  the 
stress-dependent  macrodislocation  velocity.  The  first  term  on  the  rhs.  of  Eq.  (5-18)  is  the  elastic 
strain  rate,  and  the  second  term  is  the  plastic  strain  rate  due  to  granular  slip  on  the  slip  plane. 

Next  we  operate  on  Eq.  (5-16)  with  <M)y  and  on  Eq.  (5-17)  with  c Vdt,  and  combine  with 
Eq.  (5-18)  to  obtain: 


1  32r_  1  d2r  d(NdbBvd) 
pdy2~  G  dt2  +  dt 


The  FRAGBED  2  comminution  and  dislocation  nucleation  models  yield 


l_dN±  =  _2dB  1  d(Ndvd) 
Nd  dt  B  dt  Nd  dx 


(5-20) 
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Combining  Eqs.  (5-19)  and  (5-20)  yields 


d2r  Gdh 
dt 2  p  dy2 


+A—+F  =  0 
dt 


(5-21) 


where  A  =  bBGNd  dvd/dt 
and 


F=bBGNdvd 


1  dB 
B  ddt 


i  SOW 

Nd  dx 


(5-22) 


We  first  note  that  A  is  always  positive,  because  b,  B,  G,  Nd,  and  dvd/dz  are  positive.  The 
reason  that  dvd/di  is  positive  is  that  the  model  requires  higher  shear  stresses  to  increase  the 
dislocation  velocity,  i.e.,  the  system  is  dissipative. 

The  FRAGBED  2  comminution  model  makes  dB/dt  a  function  of  the  stresses.  Thus  F  is 
a  function  of  the  stresses,  Nd,  and  the  macrodislocation  flux  in  the  x  direction,  d(NdVd  )/dx . 

We  next  linearize  Eq.  (5-21)  by  holding  all  the  variables  constant  at  their  values  at  a 
given  instant  (labeled  by  the  subscript  0)  and  then  replacing  tby  t0  +  A .  Then  Eq.  (5-21) 
becomes 


d2A 
dt 2 


Gd2  A 
pdy2 


+^+F°~° 


(5-23) 


Next  we  define 


D  =  A  +  ^ 


(5-24) 
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Combining  Eqs.  (5-23)  and  (5-24)  yields 


d2D  Gd2D  .  3D. 

W~J 


(5-25) 


Eq.  (5-25)  is  of  the  form  shown  by  Whitham  to  have  solutions  that  are  both  unique  and  stable. 
Stability  means  that  D  decreases  with  time.  Taking  the  time  derivatives  of  Eq.  (5-24)  yields: 


a d 

dt 


=  dA  +  fi 

dt  Aq 


<0 


(5-26) 


F0  can  be  positive  or  negative,  depending  on  the  comminution  rate  and  macrodislocation 
flux.  If  the  comminution  rate  is  zero,  for  example,  F0  will  be  negative  if  the  macrodislocation 
density  is  decreasing  because  more  macrodislocations  are  flowing  out  of  the  mesomechanical 
relevant  volume  element  (RVE)  than  into  it.  In  any  event,  if  F0  is  positive,  dA/dt  will  be 
negative,  and  must  be  greater  in  absolute  magnitude  than  F0fk0 .  However,  if  F0  is  negative, 

3 A /dt  can  be  positive,  but  must  have  an  absolute  value  less  than  that  of  F0fk0 . 

In  summary,  the  FRAGBED  2  model  should  satisfy  well-posedness  and  thermodynamic 
respectability  requirements  because  the  shear  wave  equations  associated  with  each  slip  plane 
satisfy  Whitham ’s  criteria.  Because  the  same  analysis  applies  to  each  slip  plane,  the  uniqueness 
and  stability  requirements  should  also  be  satisfied  for  an  RVE  containing  many  such  planes.  The 
physical  reason  for  thermodynamic  respectability  is  that  the  granule  comminution  and  flow 
process  is  dissipative. 

Internal  consistency  is  defined  as  the  requirement  that  the  mesomechanical  model  does 
not  predict  large  variations  in  calculated  variables  over  the  dimension  of  the  RVE.  Compliance 
with  both  this  requirement  and  that  of  mesh  size  insensitivity  must  be  demonstrated  for  each 
specific  application. 
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SECTION  6 

EXPERIMENTS  FOR  CALIBRATION  AND  VALIDATION  OF  FRAGBED2 


CAVITY  EXPANSION  EXPERIMENTS 

Calibration  of  the  comminution  parameters  in  FRAGBED2  requires  an  experiment  that 
applies  a  known,  relevant  stress  and  strain  history  to  a  specimen,  then  allows  the  specimen  to  be 
recovered  to  quantify  and  correlate  the  damage.  Impact  recovery  experiments  in  which  a  thin 
stelliform  flyer  impacts  a  ceramic  specimen  to  produce  a  single  well-characterized  stress  pulse 
would  be  useful  for  calibrating  the  comminution  model.  Such  experiments  have  been  performed 
on  aluminas  similar  to  our  AD-995  [20],  but  only  at  stress  levels  insufficient  to  cause  significant 
comminution. 

In  work  reported  elsewhere  [21,  22],  experiments  were  performed  by  detonating  an 
explosive  charge  within  a  spherical  cavity  machined  in  a  ceramic  specimen  confined  in  an 
impedance-matching  bronze  container.  The  resulting  loading  of  the  material  (stress,  strain,  and 
strain  rate)  is  well  characterized,  nominally  one-dimensional,  and  is  similar  to  the  conditions 
occurring  near  the  tip  of  a  rod  projectile  during  a  penetration  event.  The  test  is  instrumented  with 
electromagnetic  particle  velocity  history  gages  located  at  several  radii  from  the  charge.  Posttest 
microscopic  examination  of  the  recovered  samples  has  revealed  an  overall  damage  less  severe 
than  in  typical  penetration  experiments,  and  has  allowed  a  more  detailed  description  of  the 
fracture  and  comminution  processes. 

The  spherical  cavity  expansion  experiments  [21,22]  provide  an  order-of-magnitude 
estimate  of  fragment  size  as  a  function  of  time  and  stress  level.  In  experiments  on  AD-995, 
stresses  near  the  cavity  wall  were  estimated  to  be  around  8  GPa,  pulse  durations  were  on  the 
order  of  1  ps,  and  fragment  sizes  were  on  the  order  of  10  pm.  Corresponding  normalized  values 
are  t  =  275 ,  and  B=  1.  Comparing  these  values  with  the  results  shown  in  Figure  5,  we  see  that 
B  =  1  at  f=275  corresponds  to  a  value  of  m  greater  than  1.15. 

The  FRAGBED2  model  in  L2D  does  a  reasonably  good  job  of  matching  particle  velocity 
histories  measured  in  spherical  cavity  expansion  tests  on  Coors  AD-995  alumina  (see  Figure  6). 
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VELOCITY  (m/s)  VELOCITY  (m/s) 


Figure  6.  Comparison  of  measured  radial  particle  velocity  at  different  radii 
with  those  computed  by  FRAGBED2  in  L2D.  (continued) 
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VELOCITY  (m/s)  VELOCITY  (m/s) 


Figure  6.  Comparison  of  measured  radial  particle  velocity  at  different  radii 
with  those  computed  by  FRAGBED2  in  L2D.  (concluded) 
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The  results  are  sensitive  to  the  initial  fragment  size  distribution  B0  computed  using  Eq.  3-6.  We 
found  that  the  results  were  also  sensitive  to  the  choice  of  limiting  crack  speed  (or  breaking 
frequency),  CR  (see  Eq.  4-13).  CR  is  essentially  a  parameter  that  scales  time.  We  found  that  if  B0 
was  too  small  or  Cr  too  large,  comminution  proceeded  too  quickly,  the  response  was  too  non¬ 
elastic,  and  the  velocity  waveforms  did  not  display  the  negative  rebound  phases  that  are  in  the 
measured  histories.  By  our  choice  of  m,  B0,  and  CR,  we  could  compute  final  fragment  sizes  that 
matched  the  experiments.  In  particular,  the  variable-size  initial  fragment  bed  approach  outlined 
inEq.  3-1  to  3-6  was  used,  with  Kc  =4MPaVm,  P0  =15.5GPa,  and  a  =  0.01  m.  Those 
parameter  values,  along  with  a  uniform  initial  porosity,  <j)0  =  5.85  x  10-6 ,  m  =  1.5,  and  Cr  =  496 
m/s  (see  Eq.  4-10  through  4-15),  gave  the  result  that  nearly  all  of  the  comminution  occurred 
within  about  a  millimeter  of  the  charge,  as  was  observed  in  the  experiments.  In  the  thin 
comminuted  zone,  the  porosity  rose  rapidly  to  3%  within  6  ps,  by  the  influx  of  dislocations  from 
the  wall  of  the  spherical  cavity.  Experiments  with  stronger  charges  and  therefore  greater 
comminution  would  be  useful  for  calibrating  the  comminution  model. 

INVESTIGATION  OF  SHEAR  LOCALIZATION  IN  CERAMICS 
USING  A  MODIFIED  THICK-WALLED  CYLINDER  METHOD 

Note:  The  thick-walled  cylinder  experiments  were  designed  with  the  help  of 

Professor  Vitali  F.  Nesterenko,  University  of  California,  San  Diego.  He  will  be  a  co-author  of 

the  journal  manuscript  resulting  from  the  following  section. 

Compared  to  steels,  ceramic  armors  are  attractive  from  the  standpoint  of  weight.  To 
design  ceramic  armors  efficiently,  good  computational  models  are  needed  that  include  the 
relevant  material  physics.  Shockey  et  al.  [1]  describe  the  relevant  physics  to  include  the  fracture 
of  the  ceramic,  the  comminution  of  the  fractured  ceramic  into  a  fragmented  bed  (the  Mescall 
zone),  and  the  subsequent  flow  of  the  fragmented  bed  out  of  the  penetrator’s  path.  Curran  et  al. 
developed  a  computational  model,  FRAGBED  [2,21],  that  captures  these  physics.  Other  models 
have  also  been  developed  for  ceramic  behavior  under  armor  penetration  conditions  [23-25].  To 
calibrate  any  of  these  models,  tests  are  required  that  exercise  the  relevant  physics  in  relative 
isolation.  In  particular,  a  test  is  needed  that  subjects  ceramic  to  intense,  high  rate,  compressive 
flows  like  those  in  the  Mescall  zone,  which  is  the  zone  immediately  in  front  of  an  advancing 
penetrator.  Nesterenko  et  al.  developed  the  thick-walled  cylinder  method  (TWCM)  [26,27]  to 
meet  these  requirements. 
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Nesterenko  et  al.  applied  the  TWCM  to  study  the  large  deformation,  high  strain  rate 
behavior  of  metals  [28],  reactive  mixtures  [29],  ceramic  powders  [30],  and  intact  ceramics  [31]. 
In  the  TWCM,  depicted  schematically  in  Figure  7,  a  copper  tube  is  inserted  inside  a  thick-walled 
tube  of  specimen  material.  This  assembly  is  in  turn  inserted  inside  a  copper  outer  tube.  The 
outer  tube  is  surrounded  by  explosive  which,  when  detonated,  compresses  the  assembly  radially 
inward.  By  sectioning  the  recovered  assembly  and  measuring  the  radial  displacement,  the 
circumferential  strain  in  the  ceramic  can  be  quantified  and  the  amount  of  fragmentation  and  shear 
banding  related  to  the  strain.  The  TWCM  is  attractive  because  it  allows  specimens  to  be 
subjected  to  intense  rapid  loading  yet  be  recovered. 

In  the  original  application  of  the  TWCM  to  ceramics  [30,31],  the  tube  assembly  was 
subjected  to  two  explosions.  For  the  first  explosion,  a  solid  copper  rod  rather  than  a  tube  was 
inserted  in  the  inner  bore  of  the  specimen.  The  first  explosion  fragmented  and  compacted  the 
specimen  into  an  in  situ  fragment  bed  but  did  not  strain  it  much.  Before  the  second  explosion, 
the  solid  copper  rod  was  drilled  out  to  form  a  tube.  The  second  explosion  compressed  the 
assembly  inward  until  the  drilled  hole  collapsed  completely.  The  amount  of  strain  was  controlled 
by  varying  the  diameter  of  the  drilled  hole.  Clearly,  the  need  for  two  explosions  complicates  the 
TWCM  when  applied  to  ceramics. 

We  realized,  however,  that  the  need  for  two  explosions  can  be  eliminated.  As  shown  in 
Figure  8,  by  tapering  the  bore  of  the  inner  copper  tube  or  by  inserting  a  tapered  pin  into  a  straight 
bore,  the  diameter  of  the  bore  (and  therefore  the  allowed  range  of  strain)  can  be  varied  along  the 
length  of  the  assembly  from  near  zero  at  the  solid  end  to  as  much  as  needed  at  the  open  end.  The 
tapered  pin  scheme  was  chosen  for  implementation  because  of  the  difficulty  associated  with 
machining  a  tapered  bore. 

In  this  section  of  the  report,  we  (1)  describe  the  tapered  TWCM,  (2)  show  its  application 
to  alumina  ceramic  to  illustrate  high  rate  flow  phenomenology  and  the  transition  to  localized 
deformation,  and  (3)  present  data  for  validating  the  FRAGBED2  ceramic  comminution  and 
fragment  flow  model.  Details  of  the  tapered  TWCM  are  presented  so  that  others  can  model  the 
experiment.  Results  for  Coors  AD-998  alumina  are  shown.  Among  the  results,  particular 
emphasis  is  placed  on  shear  band  phenomenology  and  implications  for  the  fragment  flow  picture 
of  ceramic  deformation.  Also  included  are  the  results  of  simulating  the  experiment  using  a  two- 
dimensional  hydrocode. 
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Figure  7.  Original  Nesterenko  thick-walled  cylinder  method. 
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Figure  8.  Modified  thick-walled  cylinder  method.  (Continued) 


PETN  (0.67  g/cm  3) 


Figure  8.  Modified  thick-walled  cylinder  method.  (Concluded) 


Extruded  AD-998  alumina  ceramic  tubes  1 1.1 -mm  outside  diameter,  7.9-mm  inside 
diameter,  and  44.5-mm  long  were  obtained  from  Coors  Ceramics,  Golden,  CO.  The  alumina  was 
reported  to  be  99.8%  pure  with  porosity  around  1%.  The  porosity  is  estimated  from  the  ratio  of 
the  theoretical  density  of  a-alumina,  3.97  g/cm3,  to  that  of  extruded  AD-998,  3.92  g/cm3.  As 
depicted  in  Figure  8,  the  tubes  were  fitted  into  and  onto  copper  tubes  with  about  50  pm  of 
clearance.  The  assemblies  were  completed  by  inserting  copper  rings  below  the  specimens  and 
inserting  the  brass  pin  through  the  conical  cap  and  into  the  inner  tube.  Finally,  the  momentum 
trap  was  bonded  to  the  bottom  of  the  assembly  and  the  brass  pin  was  bonded  to  the  cap  with 
cyanoacrylate  adhesive.  The  assembly  was  then  placed  inside  a  100-mm  length  of  PVC  pipe  and 
topped  with  a  foam  attenuator.  About  130  g  of  PETN  explosive  was  packed  around  the  assembly 
and  covered  with  a  disk  of  sheet  explosive.  The  explosion  was  initiated  from  the  center  of  the 
disk.  The  recovered  assemblies  were  sectioned,  polished,  and  micrographed. 

Transverse  optical  photographs  of  sections  at  10-mm  intervals  measured  from  the  most 
deformed  end  of  the  ceramic  are  shown  in  Figure  9,  and  a  longitudinal  section  is  shown  in 
Figure  10.  Average  radial  displacements,  u,  finite  Lagrangian  hoop  strains,  £ee,  and  radial  strains, 
fir,  are  computed  for  each  section  from  the  radial  displacement,  uq,  of  the  inner  wall  of  the  inner 
copper  tube  from  its  initial  position,  Rq  =  2.38  mm,  to  its  final  position  against  the  tapered  pin, 
assuming  incompressible,  plane  deformation. 


u  =R--^R2  -2/^,m0  +  Mq 
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R 2  —  27?qWq  +  Uq 
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The  strains  are  plotted  in  Figure  11. 

As  revealed  in  the  Figure  9  cross  sections,  a  significant  portion  of  the  total  deformation  is 
accommodated  by  displacement  jumps  across  shear  bands.  Figure  12  shows  that  if  all  the 
deformation  were  accommodated  in  one  logarithmic  spiral  shear  band,  then  an  approximate 
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Figure  9.  TWCM  cross  sections  from  different  locations  along  the  taper,  with  average 
hoop  strains  shown,  (continued) 


Figure  9.  TWCM  cross  sections  from  different  locations  along  the  taper,  with  average 
hoop  strains  shown,  (continued) 
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Figure  9.  TWCM  cross  sections  from  different  locations  along  the  taper,  with  average 
hoop  strains  shown,  (continued) 
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Figure  9.  TWCM  cross  sections  from  different  locations  along  the  taper,  with  average 
hoop  strains  shown,  (concluded) 
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Figure  10.  Longitudinal  TWCM  cross  section. 
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HOOP  STRAIN  RADIAL  STRAIN 


LAGRANGIAN  RADIAL  POSITION  (mm) 


Figure  1 1 .  Finite  Lagrangian  strains  in  the  ceramic  for  the  cross  sections 
shown  in  Figures  9b,  c,  and  d,  assuming  incompressible,  plane 
deformation.  Inside,  outside,  and  mean  tube  radii  are  indicated. 
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relationship  between  the  circumferential  strain  and  the  displacement  jump,  S,  across  the  shear 
band  is 


AS  _  6  _  _  u 


(6-4) 


where  S  is  the  circumference  and  u  is  the  average  radial  displacement.  Thus,  when  all  of  the 
deformation  is  accommodated  by  shear  bands,  the  relationship  between  the  average  radial 
displacement  and  the  displacement  jump  across  each  of  N  bands  is 


(6-5) 


Figure  13  plots  the  displacement  and  displacement  jump  for  the  sections  in  Figure  9,  assuming 
incompressible  plane  deformation  and  assuming  that  all  of  the  deformation  is  accommodated  in 
N  =8  identical  bands.  Displacement  jumps  at  the  section  with  the  largest  strain.  Figure  9d,  are 
computed  to  be  around  700  pm. 

SEM  photographs  of  the  same  sections  shown  in  Figure  9  were  taken  of  individual  shear 
bands,  cracks,  and  relatively  undeformed  regions.  From  the  SEM  micrographs  of  the  shear 
bands,  shear  band  displacement  jumps  can  be  measured.  In  some  cases,  the  size  of  the  jump  is 
different  on  the  inside  and  the  outside  of  the  tube,  implying  that  the  material  on  one  side  of  the 
band  is  stretched.  This  is  particularly  evident  for  the  band  on  the  left-hand  side  in  Figure  14,  in 
which  there  is  virtually  no  step  on  the  outside  diameter  and  a  large  step  on  the  inside  diameter 
and  many  open  cracks  on  the  rhs  of  the  band.  Figure  15  is  a  magnified  view  of  the  cracked 
region. 

The  shear  bands  vary  in  character  independently  of  where  they  are  around  its 
circumference  or  along  the  length  of  the  tube,  i.e.,  independently  of  the  amount  of  strain.  Some 
bands  are  very  sharp  and  crack-like,  e.g.,  Figure  16;  some  are  more  broad  and  diffuse,  e.g., 
Figure  17;  some  are  difficult  to  discern  at  higher  magnification,  e.g.,  Figure  18;  some  contain 
debris  at  their  ends,  e.g.,  Figure  19;  and  some  do  not,  e.g.,  Figure  20. 
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Figure  14.  Pair  of  shear  bands  with  dilated  region  on  the  right-hand  side  of  the  left-hand  band. , 
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Figure  15.  Magnified  view  of  the  dilated  region  in  Figure  6. 
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Figure  16.  Sharp,  highly  localized  shear  bands. 
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Figure  18.  High  magnification  view  of  shear  band  at  e00  =  0.04.  Band  is  located  between  the  dashed  lines. 
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Figure  19.  Debris  accumulated  at  OD  end  of  the  shear  band  at  %  «  0.12. 
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Shear  band  quantity  is  relatively  independent  of  the  amount  of  deformation.  Although  the 
section  taken  where  the  core  was  solid  (Figure  9a)  has  no  bands,  about  5  cracks  are  oriented 
roughly  radially  and  spaced  non-uniformly.  About  10  bands  are  at  the  lowest  level  of 
compression  (Figure  9b),  with  about  4  of  those  having  significant  shear  steps  on  the  ID  or  OD. 
About  8  bands  are  at  the  next  level  of  compression  (Figure  9c),  all  having  significant  shear. 
Finally,  about  9  bands  are  at  the  highest  level  of  compression  (Figure  9d)  with  even  more 
significant  shear  in  each.  Thus,  the  number  of  bands  appears  to  be  roughly  constant  at 
8-10  spaced  non-uniformly  in  the  circumferential  direction  independently  of  the  amount  of 
inward  compression. 

Two-dimensional,  free-Lagrange,  hydrocode  simulations  of  the  experiment  were 
performed  using  the  FRAGBED2  implementation  in  the  L2D  hydrocode  [13]  to  estimate  the 
stress  and  strain  rate  conditions  in  the  ceramic.  The  model  included  the  explosive  detonation. 
Two  sets  of  simulations  were  performed.  In  one  set,  cylindrical  axisymmetry  was  assumed,  that 
is,  the  problem  coordinates  were  the  axial  and  radial  directions.  Figure  21  shows  the  mesh  and 
Table  3  lists  input  parameters.  The  copper  was  modeled  as  Von  Mises  elastic-plastic.  The 
axisymmetric  model  cannot  capture  the  shear  banding  behavior  and  the  simulation  therefore 
provides  only  qualitative  estimates  of  the  average  conditions  in  the  ceramic. 

Table  3 

FRAGBED2  TWCM  AND  BALLISTIC  SIMULATION  PARAMETERS 


Parameter 

Symbol 

Value 

Units 

Fracture  Toughness 

Klc 

4 

MPaVm 

Burgers  Vector 

b 

0.01 

Initial  Block  Size" 

B0 

3 

pm 

Initial  Dislocation  Density 

No 0 

1.14x10'° 

nf2 

Comminution  Rate  Saturation  Width 

m 

1.5 

Relative  Crack  Size 

n 

0.01 

Max  Crack  Speed 

496 

m/s 

Friction  Coefficient 

V 

0.3 

±20%  random  variation  introduced  in  plane  strain  simulation. 


In  the  second  set  of  simulations,  plane  strain  deformation  was  assumed  and  the  problem 
coordinates  were  the  radial  and  circumferential  directions  (see  Figure  22  and  Table  3).  The 
simulation  corresponded  to  modeling  the  behavior  of  a  single  axial  location,  e.g..  Figure  9d.  As 
in  the  axisymmetric  simulations,  the  copper  was  modeled  as  elastic-plastic,  and  the  alumina  was 
modeled  using  the  FRAGBED2  model.  The  FRAGBED2  model  uses  an  analogy  to  dislocation 
theory  to  describe  the  deformation  of  the  ceramic.  Fragments  are  imagined  to  slide  on  discrete 
planes  at  speeds  controlled  by  the  local  stresses.  Fragments  can  also  reduce  in  size  according  to  a 
comminution  rate  law.  In  an  attempt  to  precipitate  localization,  a  random  ±  20%  variation  in 
initial  fragment  size  and  dislocation  density  was  introduced  in  the  axisymmetric  simulations. 

The  results  from  the  axisymmetric  model  show  that  the  radial  stress  pulse  is  compressive 
and  persists  for  about  10  ps  and  peaks  around  4-8  Gpa;  the  circumferential  stress  pulse  plateaus 
at  about  the  input  Hugoniot  elastic  limit  for  the  ceramic,  6  Gpa,  and  has  a  duration  of  up  to 
15  ps;  and  the  axial  stress  is  compressive,  about  2  Gpa  in  magnitude,  and  persists  for  about 
25  ps.  The  radial  velocity  is  around  1.5  m/s  inward.  Radial  strains  are  overall  positive  and  vary 
from  0%  to  20%  along  the  taper  section,  and  the  radial  strain  rate  is  about  2  x  104  s'1.  Hoop 
strains  are  compressive  and  vary  up  to  -15%  at  similar  rates.  Axial  strains  are  at  most  5%. 

Figure  23  shows  strain  paths  in  the  r-0  plane. 

The  fact  that  the  radial  stress  pulse  is  5  ps  shorter  than  the  circumferential  stress  pulse 
and  10  ps  shorter  than  the  axial  stress  pulse  may  explain  some  of  the  propensity  for  shear 
banding  and  the  radially  stretched  region  on  the  right  side  of  the  band  in  Figures  14  and  15.  The 
circumferential  and  axial  stresses  remain  high  after  the  radial  stress  falls  off,  which  places  the 
ceramic  in  a  state  of  biaxial  compression.  The  reduced  stress  in  the  radial  direction  implies  that 
the  ceramic  might  be  able  to  expand  in  the  radial  direction  and  that  faulted  blocks  might  be  able 
to  slide  off  each  other  in  the  radial  direction. 

The  deformation  computed  in  the  plane  strain,  FRAGBED2  model  shows  only  slight 
tendency  to  localize  into  logarithmic  spiral  shear  bands.  As  shown  in  Figures  24  and  25,  there 
are  indications  of  localized  comminution  and  localized  deformation,  but  the  localization  is  of  a 
different  form  and  more  diffuse  than  observed  in  the  experiments,  and  might  be  related  to  the  fact 
that  fragment  dislocations  are  restricted  to  move  on  a  set  of  discrete  planes,  rather  than  the 
emergence  of  a  true  softening  instability.  It  might  also  be  related  to  the  fact  that  in  the  plane 
strain  simulation  there  is  a  strong  radial  shock,  whereas  in  the  TWCM  test  the  wave  sweeps 
down  the  length  of  the  ceramic  tube.  In  the  simulation,  the  radial  motion  rebounds  strongly  after 
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Figure  22.  L2D  mesh  for  plane  strain  FRAGBED2  simulations  of  the  TWCM. 

Scale  in  centimeters.  Numbers  with  asterisks  indicate  time  history  data  plot  points. 
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the  inner  copper  tube  contacts  the  central  pin,  but  there  is  no  evidence  of  rebound  in  the  TWCM 
assemblies  sectioned  posttest  (see  Figures  9  and  10).  The  sweeping  motion  combined  with  the 
shear  localization  makes  the  TWCM  fully  three-dimensional.  Thus,  plane  strain  simulations  are 
not  completely  faithful  representations  of  the  true  deformation.  Perhaps  the  sweeping  motion 
helps  trigger  localization.  On  the  other  hand,  the  experimental  rock  mechanics  literature 
contains  instances  of  quasistatic  and  dynamic  plane  strain  radial  compression  of  holes  leading  to 
spiral  shear  bands. 

The  transverse  sections  show  clearly  that  shear  localization  or  faulting  is  a  significant 
mechanism  by  which  the  hoop  compression  is  accommodated,  even  at  the  smallest  levels  of 
hoop  compression  (Fig.  9b).  However,  the  displacement  jumps  in  Figure  9  are  significantly 
smaller  than  those  estimated  as  if  all  the  deformation  were  by  shear  banding  (Figure  13),  which 
implies  that,  indeed,  not  all  of  the  deformation  is  accommodated  by  the  bands.  Another 
significant  mechanism  must  be  present  for  accommodating  hoop  compression.  Given  the 
fragmentation  visible  in  Figure  9  and  10,  and  given  the  low  amount  of  atomic  dislocation 
plasticity  observed  in  other  explosively  loaded  alumina  samples  [32],  shifting,  rearrangement, 
and  breaking  of  fragments  appears  to  be  a  likely  mechanism. 

The  transverse  section  taken  where  the  core  was  solid  and  where  there  was  therefore  only 
minimal  straining  shows  radial  cracks  (Figure  9a)  that  could  immediately  lead  to  the  shear 
banding  mechanism  with  further  compression.  For  instance,  a  detail  of  one  of  these  cracks. 

Figure  26,  shows  that  it  is  curved  and  that  further  compression  could  be  accommodated  by 
sliding  along  the  crack  faces. 

Chen  and  Ravichandran  [33]  performed  split-Hopkinson  bar  compression  tests  on  glass- 
ceramic  cylinders  surrounded  by  thick  metal  jackets  intended  to  provide  radial  confining  stress. 

All  of  their  confined  cylinders  failed  by  faulting,  and  the  faulting  appeared  to  begin  at  the  earliest 
stages  of  the  deformation.  As  was  the  case  for  TWCM,  there  seemed  to  be  little  or  no 
homogeneous  inelastic  deformation  prior  to  the  onset  of  faulting.  In  contrast  with  TWCM, 
however,  the  amount  of  debris  in  the  fault  slip  region  seemed  to  correlate  with  the  amount  of 
slip. 

Nesterenko  et  al.  applied  the  untapered  TWCM  to  compacted  ceramic  powders  [30]. 
Compared  with  the  results  obtained  here,  the  powder  tests  produce  many  more  bands  and  the 
bands  are  more  uniform  in  displacement  jump  and  spacing.  The  same  is  true  for  Nesterenko’s  test 
on  ceramic  subjected  to  two  explosions.  In  either  case,  more  than  50  bands  were  produced  at  an 
average  circumferential  strain  around  -0.2.  As  in  the  work  reported  here,  Nesterenko  et  al.  [30] 
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Figure  26.  Curved  crack  in  section  with  solid  core,  f00 
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found  that  significant  portions  of  the  total  deformation  were  accommodated  by  mechanisms  other 
than  shear  banding. 

The  flowing  fragment  bed  picture  of  ceramic  deformation  differs  from  the  shear 
localization  picture  observed  in  the  TWCM  tests.  Which  picture  applies  to  ceramic  armor 
penetration  is  an  important  issue.  Apparently,  both  pictures  apply,  perhaps  at  different  stages  in 
the  penetration  process.  This  issue  has  strong  implications  for  armor  design,  modeling,  and 
choice  of  experiments  for  characterizing  penetration  resistance. 

The  preponderance  of  shear  bands  and  the  fact  that  even  the  smallest  levels  of  radial 
compression  result  in  cracks  that  could  become  bands  suggest  that,  under  the  geometry,  stresses, 
strain,  and  loading  rate  conditions  under  consideration,  shear  banding  is  a  significant  mechanism 
by  which  the  ceramic  deforms.  However,  fragment  flow  also  appears  to  be  a  significant 
mechanism. 

BALLISTIC  TESTS 

Subramanian  and  Bless  [34]  have  performed  a  series  of  ballistic  experiments  in  which 
they  shot  confined  AD-995  alumina  cylinders  against  stationary  L/D=20  tungsten  rod 
penetrators  in  the  reverse  ballistics  mode.  Penetration  and  rod  erosion  histories  were  recorded 
using  flash  x-rays.  We  simulated  their  tests  at  1.5  and  3.5  km/s  using  FRAGBED2  in  L2D,  with 
the  same  set  of  parameters  for  both  simulations.  Parameter  values  are  shown  in  Table  3. 

Figure  27  shows  the  computational  mesh  at  a  stage  part  way  through  a  penetration.  Figure  28 
shows  the  computed  penetrator  tip  and  tail  trajectories  for  the  two  simulations,  along  with 
average  steady-state  velocities  reported  by  Subramanian  and  Bless  [34].  At  1.5  km/s,  the 
computed  penetration  rate  is  well  within  the  scatter  band  of  the  measured  rates.  At  3.5  km/s, 
computed  rates  are  barely  within  the  low  side  of  the  scatter  band,  about  15%  below  the  mean  of 
the  measured  rates. 

The  computed  results  are  sensitive  to  initial  fragment  size  and  to  the  initial  porosity. 
Subramanian  and  Bless  claim  an  initial  porosity  for  their  AD-995  of  around  2%.  Our 
microscopic  observations  of  AD-995  indicate  a  porosity  well  below  5%.  Based  on  the  ratio  of 
the  manufacturer’s  reported  density  for  AD-995, 3.89  g/cm3,  to  that  of  a-alumina,  3.97  g/cm3, 
we  also  estimate  the  porosity  to  be  around  2%.  In  any  case,  we  found  that  FRAGBED2,  with  an 
initial  porosity  around  2%,  gave  the  best  match  to  the  ballistic  experiments  of  Subramanian  and 
Bless. 
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Figure  27.  FRAGBED2/L2D  simulation  of  tungsten  long  rod  penetration  in  AD-995 
alumina.  Scale  is  in  centimeters. 
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SECTION  7 

DISCUSSION  AND  CONCLUSIONS 


In  this  section,  we  list  the  FRAGBED2-specific  input  parameters  and  discuss  the 
importance  of  each  in  determining  key  aspects  of  penetration  or  cratering  behavior.  FRAGBED2 
requires  the  following  FRAGBED2-specific  input  parameters.  (Non-FRAGBED2-specific 
parameters  include  the  solid  density  and  the  solid  moduli.) 


Table  4 

FRAGBED2  PARAMETERS 


Parameter 

Description 

Comments 

No.  and  orientation  of  planes 

Arbitrary;  can  be  tailored  to 
each  problem 

All  non-elastic  slip  occurs  on 
these  planes. 

B0( cm) 

Characteristic  fragment  granule 
size  at  onset  of  granular  flow 
and  comminution 

This  parameter  is  determined 
from  fracture  calculations  or 
algorithms. 

A/oo  (no./cm2) 

Initial  dislocation  density 
(number/unit  area) 

Depends  on  boundary  and  initial 
conditions 

b  \ 

Dislocation  jog  coefficient 

Related  to  granule  morphology 

Coefficient  of  friction  between 
granules 

Normally  based  on  handbook 
values 

h,  m 

Comminution  parameters 

Chosen  by  correlation  with 
cavity  expansion  experiments 

9 

Granule  geometry  factor 

Normally  set  to  unity 

K1D(dyne  cm'3'2) 

Solid  dynamic  fracture  initiation 
toughness 

Normally  based  on  handbook 
value 

mi,  m2 

Local  stress  factor  coefficients 

Related  to  granule  morphology 

All  ten  parameters  in  the  above  table  can  be  in  principle  assigned  different  values  on  each 
plane  of  this  multiplane  model.  (This  is  in  fact  an  important  capability  for  cases  in  which  pre- 
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existing  interfaces  or  planes  of  weakness  exist  in  the  material.)  The  total  number  of  potential 
input  parameters  is  thus  equal  to  10  times  the  number  of  planes  used  in  the  model.  For  example, 
in  the  L2D  applications  discussed  above  we  used  1 3  planes,  so  the  number  of  potential 
parameters  in  that  case  is  130. 

However,  in  initially  isotropic  cases,  we  use  the  same  parameters  on  all  planes.  Thus,  in 
many  if  not  most  armor  ceramic  applications  (such  as  the  examples  shown  in  this  report),  the 
number  of  FRAGBED-specific,  input  material,  property  parameters  is  10. 

The  effects  of  varying  input  parameters  in  FRAGBED2  are  difficult  to  predict  a  priori 
because  of  the  many  interactions  between  different  parts  of  the  model.  In  a  parametric  study 
leading  to  the  correlation  between  computations  and  experimental  data  discussed  above,  we 
varied  all  the  above  input  parameters  except  g  (assumed  to  be  unity)  and  K,d  (set  to 
4  MPa^/n~  y  The  results  showed  that,  as  hoped,  they  depended  on  the  input  parameters  in  a 
stable  way,  i.e.,  modest  changes  in  the  input  cause  only  modest  changes  in  the  output. 

In  conclusion, 

•  The  FRAGBED2  model  treats  the  physics  relevant  to  ceramic  armor 
penetration:  fracture,  comminution,  and  flow  of  the  fragmented  bed. 

•  FRAGBED2  models  reasonably  well  laboratory  experiments  in  which  alumina 
armor  was  impacted  by  tungsten  long  rods. 

•  FRAGBED2  correctly  models  particle  velocity  histories  obtained  in  spherical 
cavity  expansion  experiments. 

•  FRAGBED2  predictions  are  sensitive  to  readily  identifiable  and  measurable 
material  parameters,  such  as  initial  fragment  size,  fracture  toughness,  and 
friction  coefficient. 

•  We  have  not  found  a  set  of  FRAGBED2  parameters  that  would  allow  it  to 
predict  the  shear  bands  observed  in  TWCM  tests  performed  on  alumina. 

•  The  deformation  in  TWCM  tests  on  alumina  is  accompanied  by  both  shear 
bands  and  by  fragment  flow. 
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APPENDIX:  NUMERICAL  DEFINITIONS  AND  ROUTINES 


DEFINITIONS 

Block  size  =  initial  block  size  in  centimeters.  There  is  only  one  initial  block  size,  common  for  all 
planes.  No  default. 

Burger  =  burgerfs  vector,  non-dimensional.  No  default. 

Coheslip  =  initial  cohesion  on  all  slip  planes  in  MBar.  Default  zero. 

Crayleig  =  Rayleigh  wave  velocity  in  cm/|Lis.  No  default. 

2 

Dislocde  =  initial  dislocation  density  on  all  planes  in  cm'  .  No  default. 

Elastic:  elastic  only  calculation.  No  slip  on  planes.  For  testing  purposes. 

Endsub:  end  of  FRAGBED  input. 

Eta  =  T|  in  the  comminution  model,  dimensionless.  No  default. 

Ftoughne  =  fracture  toughness  in  MBar  |  cm.  No  default. 

Gsmall  =  g  in  the  formula  for  the  plastic  strain  rate  on  a  slip  plane,  dimensionless.  No  default. 
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Initfp:  non-trivial  initialization  of  the  plastic  deformation  tensor  Fp.  Ordinarily  Fp  is  initialized 
to  unity.  This  flag  will  make  the  code  initialize  Fp  to  a  non-trivial  value  that  is  hard-wired  in 
the  code  (look  for  the  string  NON-TRIVIAL  FP).  For  testing  purposes. 

Ml  =  the  shear  stress  factor  in  the  driving  stress  formula,  dimensionless.  No  default. 

M2  =  the  normal  stress  factor  in  the  driving  stress  formula,  dimensionless.  No  default. 

Mufricti  =  the  coefficient  of  friction  on  all  planes,  dimensionless.  No  default. 

Notslip:  ordinarily  two  pairs  of  out-of-plane  planes  allow  slip  that  limits  the  out-of-plane  shear 
stresses.  This  flag  will  limit  the  planes  to  the  ones  that  the  user  specifies  in  the  x-y  plane. 
Setting  NOTSLIP  in  a  2D  problem  would  cause  all  slip  to  occur  on  planes  with  normals  in 
the  x-y  planes.  This  would  result  in  unlimited  shear  strength  in  the  two  tangential  directions. 
For  testing  purposes  and  true  3D  problems.  For  3D  simulations,  the  user  should  always  set 
the  planes  manually. 

Notransp:  turns  off  the  dislocation  transport  on  all  planes.  For  testing  purposes. 

Nslip  =  number  of  slip  planes.  In  2D,  these  are  the  planes  that  the  user  specifies.  Do  not  include 
the  default  out-of-plane  planes.  In  3D,  remember  to  set  the  notslip  flag  above.  No  default. 

Satwidth  =  the  m  parameter  in  the  saturation  and  threshold  stress  formulas  in  the  comminution 
model,  dimensionless.  No  default. 

Snormal(I)  =  Nx  Ny  Nz 

I  =  plane  number  (1  through  nslip).  No  default. 

(Nx,Ny,Nz)  is  the  slip  plane  normal.  The  code  will  normalize  the  vector.  No  default.  Always 
three  values.  In  2D,  Nz  should  be  zero. 
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ROUTINES 


SUBROUTINE  FRAGB(ffiRR,LLnjT, NCALL, DELTC,EQ,SDEV,FD,FPI,EXC,HEPS, 

PCUT,S,EXM, LABEL, RLABEL,JPOSN,LM3,NUM) 

This  is  the  main  subroutine  for  the  FRAGBED  model.  This  is  the  only  routine  that  the 
simulation  code  needs  to  interface.  There  are  essentially  three  times  when  the  simulation  code 
needs  to  call  FRAGB:  at  material  model  parameter  setup,  at  internal  state  variable  setup,  and  at 
calculation  of  the  stress  from  the  deformation.  See  NCALL  below. 

FRAGB  Arguments 

In  the  following  INTEGER  and  REAL*8,  refer  to  FORTRAN  programming  syntax.  REAL*8 
is  called  double  precision  on  most  computers. 

INTEGER  IERR:  negative  for  error.  Initially  set  to  0  by  FRAGB.  Returns  storage  used  by  the 
EXM  and  EXC  arrays. 

INTEGER  LUUT:  FORTRAN  logical  unit  number  for  the  general  output  file.  Should  be  set  at 
any  call. 

INTEGER  NCALL: 

0  for  material  model  setup.  The  material  model  parameters  must  be  supplied  via  LABEL, 
RLABEL,  etc.  for  this  call.  The  number  of  REAL*8  words  needed  for  the  material  model 
parameters  in  the  array  EXM  is  returned  in  IERR. 

1  for  initialize  internal  state  variables.  The  number  of  REAL*8  words  needed  for  the  material 
model  parameters  in  the  array  EXC  is  returned  in  IERR. 

2  for  compute  new  stress  tensor  given  new  deformation  tensor  FD.  This  is  the  critical  point 
for  this  model.  Everything  in  the  finite  difference  update  scheme  is  updated  in  this  call. 
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However,  the  dislocation  transport  requires  one  immediately  preceding  call  to  load  transport 
variables  and  one  immediately  succeeding  call  to  unload  the  updated  transport  variables. 

3  for  return  slip  plane  variables  in  the  array  EQ.  For  this  call,  EQ  must  be  an  array  of  at  least 
dimension  NSLIP*7+2.  Where  NSLEP  is  the  maximum  number  of  slip  planes  possible  in  the 
model  (not  the  particular  number  of  slip  planes  for  this  case).  This  call  can  be  used  to  extract 
information  for  contour  plots,  time  histories,  etc. 

4  for  load  transport  variables.  When  updating  the  dislocation  transport,  the  scheme  must  have 
information  about  the  internal  state  variables  for  the  cell  neighbors  to  the  cell  that  we  are 
currently  processing.  Also  needed  are  the  normal  vector  and  surface  area  of  the  interface 
between  the  current  cell  and  the  respective  neighbor.  This  information  about  the  neighbors  is 
transferred  via  the  RLABEL  and  LM3  variables  borrowed  for  this  purpose.  The  information 
is  stored  internally  in  the  subroutine  in  static  variables.  This  call  must  be  immediately 
preceding  the  compute  call  above  (NCALL=2). 

5  for  unload  transport  variables.  The  dislocation  density  is  returned  for  all  the  neighbor  cells 
that  have  been  updated.  The  dislocation  density  for  the  current  cell  itself  is  returned  as  part 
of  the  compute  call  above  (NCALL=2).  This  call  must  immediately  succeed  the  compute  call 
above  (NCALL=2). 

6  for  return  plastic  deformation  tensor  rigid  body  angle.  This  is  part  of  the  rezoning  calls 
described  below.  The  rezoning  is  by  necessity  very  code  dependent.  Also,  the  call  only 
applies  to  2D,  because  in  3D  there  are  three  angles. 

7  for  backcalculate  Fe  from  F  and  Fp  for  rezoning  purposes. 

8  for  do  the  rezoning. 

9  for  after  rezone  clean  up. 

10  for  printing  a  summary  of  the  internal  state  variables  on  the  file  connected  to  logical  unit 
LUUT.  The  printout  is  80  column  wide.  There  is  one  printout  for  every  call  to  FRAGB,  i.e., 
the  caller  must  loop  over  all  the  cells  for  which  a  printout  is  desired. 

REAL*8  DELTC:  the  current  timestep  of  the  finite  difference  scheme. 
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REAL*8  EQ:  array  whose  dimension  varies  with  the  call.  The  general  use  is  to  provide 

information  for  the  pressure  equation  of  state.  The  details  depend  on  the  model.  In  the  L2D 
code,  EQ  provides  specific  volume  and  internal  energy  and  returns  pressure,  temperature,  and 
sound  velocity.  EQ  is  also  used  for  returning  slip  plane  variables  (NCALL=3). 

REAL*8  SDEV:  array  of  dimension  five.  Transfers  stress  deviator  components.  The 
components  are  in  the  order  listed  below: 

1 :  aexx 

2.  Qe^y 

3:  aexv 

4:  aezx 

5  . 

Together  with  the  pressure  in  EQ  above  we  have  a  total  of  six  components. 

REAL* 8  FD:  array  of  dimension  3  by  3.  The  total  deformation  tensor. 

REAL*8  FPI:  array  of  dimension  3  by  3.  The  inverse  of  the  plastic  deformation  tensor. 

REAL*8  EXC:  array  of  dynamic  dimension.  The  internal  state  variables  are  stored  in  this  array. 
When  the  internal  state  variables  are  initialized  (NCALL=1),  the  total  space  in  REAL*8 
words  used  up  in  EXC  is  returned  in  IERR. 

REAL*  8  HEPS:  returns  Hillfs  equivalent  plastic  strain. 

REAL*8  PCUT:  pressure  tension  cutoff.  The  use  of  this  argument  and  S  below  depends  on  the 
elastic  model. 
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REAL* 8  S:  array.  The  dimension  depends  on  the  elastic  model.  Contains  parameters  for  the 
elastic  model.  For  L2D  S  contains  the  parameters  for  pressure  equation  of  state  plus  the 
shear  modulus  for  the  stress  deviator  model. 

REAL* 8  EXM:  array  of  dynamic  dimension.  The  parameters  for  the  slip  plane  model  are  stored 
in  this  array.  When  the  model  parameters  are  read  in  (NCALL=0),  the  total  space  in 
REAL*8  words  used  up  in  EXM  is  returned  in  IERR. 

CHARACTER*8  LABEL:  array  of  dynamic  dimension.  The  user  input  labels  for  the  material 
parameters  (NCALL=0)  are  sent  to  FRAGB  in  this  array. 

REAL* 8  RLABEL:  array  of  dynamic  dimension.  The  values  corresponding  to  the  LABELS 
above  are  sent  to  FRAGB  in  this  array. 

INTEGER  JPOSN:  array  of  dynamic  dimension.  Indicators,  indices,  and  the  like  associated  with 
LABEL  and  RLABEL  are  sent  to  FRAGB  in  this  array.  For  example,  the  following  user 
input,  “SNORMAL(2)  =  0.3”  would  result  in  LABEL  =  “SNORMAL”,  RLABEL  =  0.3,  and 
JPOSN  =  2. 

INTEGER  LM3:  the  array  position  in  LABEL,  RLABEL,  and  JPOSN  that  is  now  processed. 
Ordinarily  LM3  is  1  upon  entry  of  FRAGB.  However,  if  the  FRAGBED  input  is  part  of  a 
larger  free  format  input,  LM3  might  be  non-trivial. 

INTEGER  NUM:  the  maximum  value  for  LM3,  i.e.,  the  total  number  of  free  format  input  items. 

The  FRAGB  routine  reads  in  the  material  model  parameters,  initializes  the  internal  state 

variables,  and  does  part  of  the  stress  computations. 

The  first  step  in  the  stress  computation  is  to  call  routine  FRAGE  that  calculates  the  elastic 

stresses  given  a  certain  plastic  deformation  tensor.  Next,  routine  FRAGY  is  called  to  check  for 
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slip  on  the  slip  planes.  If  there  is  slip,  FRAGY  iterates  for  a  plastic  deformation  tensor  that 
fulfills  the  slip  condition  on  all  planes.  During  this  iteration,  FRAGE  is  called  repeatedly  from 
within  FRAGY. 

The  slip  will  limit  shear  stresses  but  will  not  limit  hydrostatic  tension.  The  routine  FRAGC  is 
called  to  check  the  normal  stress  for  tension  on  each  plane.  If  the  normal  tension  exceeds  the 
cohesion  on  the  plane,  the  pressure  is  reset  to  push  up  the  normal  stress  on  the  plane.  The 
pressure  is  not  reset  if  it  is  already  compressive.  The  slip  adjusts  the  stress  deviators,  so  there 
should  only  be  need  for  adjustment  of  the  pressure  if  the  pressure  is  tensional. 

Thereafter  the  comminution  equation  is  solved.  Finally  the  dislocation  transport  is  done.  The 
neighbor  cell  values  needed  for  the  transport  were  stored  in  arrays  local  to  FRAGB  by  a  previous 
call  to  load  the  transport  variables  (NCALL=4). 

FRAGB  also  responds  to  various  service  calls,  e.g.,  NCALL=3,  6,  10  above. 

Hillfs  equivalent  plastic  strain  is  ordinarily  defined  as  the  summation  of  the  square  root  of  the 
second  invariant  of  the  plastic  strain  increment.  Having  a  finite  deformation  tensor,  we  define 
Hillfs  equivalent  plastic  strain  instead  as  the  square  root  of  the  second  invariant  of  the  plastic 
strain  tensor. 

The  comminution  constitutes  an  ordinary  differential  equation  with  one  unknown  variable,  the 
block  size,  and  one  independent  variable,  the  driving  stress.  The  problem  with  solving  the 
equation  numerically  is  that,  for  a  certain  range  of  the  driving  stress,  the  block  size  varies  very 
rapidly.  For  other  ranges  of  the  driving  stress,  it  does  not  vary  at  all.  Because  we  know  that  the 
block  size  varies  monotonically  (it  always  decreases),  a  solution  method  with  a  variable  time  step 
is  the  best  approach. 

When  the  block  size  varies,  we  must  subiterate  on  the  comminution  equation  over  the  global 
time  step  (of  the  Lagrangian  simulation).  We  start  with  a  time  step  equal  to  the  global  time  step. 
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We  take  a  trial  increment.  If  the  block  size  decreases  by  more  than  a  factor  of  two,  we  decrease 
the  time  step  by  a  factor  of  two  and  try  again.  If  the  decrease  in  the  block  size  is  less  than  0.1%, 
we  have  a  good  increment  and  we  take  the  next  trial  increment. 

If  the  change  in  the  block  size  is  between  those  two  limits,  we  calculate  the  mean  of  the  old  and 
new  block  size,  Bmn,  and  use  this  mean  to  calculate  the  new  driving  stress  and  the  new  width  of 
the  transition.  The  increment  in  the  block  size  is  a  function  of  these  two  variables.  We  therefore 
want  the  driving  stress  and  the  transition  width  to  be  time  centered  between  the  old  and  the  new 
time  step  level.  We  then  go  back  and  increment  the  block  size  with  these  mean  values.  We 
finally  check  that  the  resulting  block  size  is  smaller  than  or  equal  to  Bmn  above.  If  it  is  not,  we 
cut  the  time  step  in  half  and  do  the  trial  increment  over  again.  This  maniac  check  is  important 
when  the  block  size  varies  very  quickly. 

If  for  a  good  increment  the  decrease  in  the  block  size  is  less  than  1%,  we  increase  the  time  step 
by  50%  on  the  next  trial  increment. 
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